Load Libraries

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
library(cowplot)
## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp
library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(ggspatial)
library(ggmap)
## ℹ Google's Terms of Service: <https://mapsplatform.google.com>
##   Stadia Maps' Terms of Service: <https://stadiamaps.com/terms-of-service/>
##   OpenStreetMap's Tile Usage Policy: <https://operations.osmfoundation.org/policies/tiles/>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
## 
## Attaching package: 'ggmap'
## 
## 
## The following object is masked from 'package:cowplot':
## 
##     theme_nothing
library(ggsn)
## Loading required package: grid
library(patchwork)
## 
## Attaching package: 'patchwork'
## 
## The following object is masked from 'package:cowplot':
## 
##     align_plots

Load data and metadata

#load wax data
wax.data <- read.csv("Data/test_bulk_density_data.csv", header=TRUE)
colnames(wax.data) <- c("Sample.ID", "DW.clean.g", "DW.wax.g", "BW.wax.g", "Temperature.C", "Notes")
wax.data$Sample.ID <- as.character(wax.data$Sample.ID)
wax.data <- na.omit(wax.data)
nrow(wax.data)==length(unique(wax.data$Sample.ID))
## [1] FALSE
duplicated_names <- wax.data[duplicated(wax.data$Sample.ID), ]
print(duplicated_names)
##      Sample.ID DW.clean.g DW.wax.g BW.wax.g Temperature.C
## 45        3685      8.470    8.965    4.824         20.09
## 121       3637      9.671   10.224    5.308         20.06
## 146       3838      6.228    6.550    3.018         20.06
## 187       3834      5.672    6.031    3.264         20.03
## 861       2050      6.196    6.566    0.515         20.01
## 941        938      8.011    8.469    4.117         20.01
## 1002       637      1.910    2.040    1.052         19.98
##                              Notes
## 45   duplicated, Pierrick has page
## 121             duplicated, page74
## 146  duplicated, Pierrick has page
## 187  duplicated, Pierrick has page
## 861           duplicated, page 101
## 941           duplicated, page 104
## 1002
species.id <- read.csv("Data/test_Haplotype_data.csv", header=TRUE)
species.id$Species <- factor(species.id$Species)
species.id$Sample.ID <- as.character(species.id$Sample.ID)

metadata <- read.csv("/Users/hputnam/MyProjects/Pocillopora-Lagoon-Abundance/Data/POC_ID_RA_list.csv", header=TRUE) %>%
  select(c("Number", "Stage", "Timing", "Site", "Transect"))

colnames(metadata) <- c("Sample.ID", "Stage", "Timing", "Site", "Transect")
metadata$Sample.ID <- as.character(metadata$Sample.ID)

#Calculate surface area (Volume enclosed) using the standard curve

#View the range of DW.clean.g
range(wax.data$DW.clean.g)
## [1]  0.023 31.154
boxplot(wax.data$DW.clean.g)

hist(wax.data$DW.clean.g)

#View the range of DW.wax.g
range(wax.data$DW.wax.g)
## [1]  0.032 32.007
boxplot(wax.data$DW.wax.g)

hist(wax.data$DW.wax.g)

wax.data <- wax.data %>%
  mutate(Vol.enc = DW.wax.g - BW.wax.g) %>%
  mutate(Bulk.Density = DW.clean.g/Vol.enc)

#check the range of bulk density values
#1.44 - 2.22 reported 10.7717/peerj.3191
range(wax.data$Bulk.Density)
## [1] -0.2958015  3.7469235
boxplot(wax.data$Bulk.Density)

hist(wax.data$Bulk.Density)

Set species colors

species_colors <- c("P. tuahiniensis" = "#D55E00",
                    "P. meandrina" = "pink",
                    "P. meandrina P. grandis" = "#0072B2",
                    "P. verrucosa" = "#E69F00",
                    "P. grandis" = "#56B4E9",
                    "P. effusa" = "#009E73",
                    "P. acuta" = "#e63946")

Plot by species and site

data <- left_join(wax.data, species.id, by= "Sample.ID")
## Warning in left_join(wax.data, species.id, by = "Sample.ID"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 4 of `x` matches multiple rows in `y`.
## ℹ Row 515 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
data <- left_join(data, metadata, by= "Sample.ID")
## Warning in left_join(data, metadata, by = "Sample.ID"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 4 of `x` matches multiple rows in `y`.
## ℹ Row 842 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
data$Species <- factor(data$Species, levels = c("P. tuahiniensis", 
                                                "P. meandrina", 
                                                "P. meandrina P. grandis", 
                                                "P. verrucosa", 
                                                "P. grandis", 
                                                "P. effusa", 
                                                "P. acuta"))

str(data)
## 'data.frame':    1009 obs. of  20 variables:
##  $ Sample.ID            : chr  "3968" "3847" "3841" "3834" ...
##  $ DW.clean.g           : num  30.62 6.59 9.7 5.73 5.73 ...
##  $ DW.wax.g             : num  31.88 6.96 10.22 6.16 6.16 ...
##  $ BW.wax.g             : num  16.16 3.35 5.03 2.89 2.89 ...
##  $ Temperature.C        : num  20.1 20.2 20.2 20.1 20.1 ...
##  $ Notes.x              : chr  "" "" "" "duplicated, Pierrick has page" ...
##  $ Vol.enc              : num  15.72 3.6 5.19 3.27 3.27 ...
##  $ Bulk.Density         : num  1.95 1.83 1.87 1.75 1.75 ...
##  $ Sanger.Plate.Name    : chr  NA "RAA_241009" "RAA_241009" "RAA_241007B" ...
##  $ Sequence.ID          : chr  NA "RA0501" "RA0492" "RA0304" ...
##  $ Plate.Well           : chr  NA "B4" "A3" "C3" ...
##  $ gDNA.Extraction.Plate: chr  NA "Plate015" "Plate015" "Plate013" ...
##  $ Haplotype            : chr  NA "Haplotype 1a" "Haplotype 1a" "Haplotype 1a" ...
##  $ Notes.y              : chr  NA "" "" "" ...
##  $ Species              : Factor w/ 7 levels "P. tuahiniensis",..: NA 3 3 3 3 3 3 2 3 3 ...
##  $ X                    : chr  NA "" "" "" ...
##  $ Stage                : chr  "" "" "" "" ...
##  $ Timing               : chr  "Post_bleach" "Pre_bleach" "Pre_bleach" "Pre_bleach" ...
##  $ Site                 : chr  "S15" "S14" "S14" "S14" ...
##  $ Transect             : chr  "T2" "T2" "T2" "T2" ...
##  - attr(*, "na.action")= 'omit' Named int [1:9] 711 712 713 714 715 716 717 718 719
##   ..- attr(*, "names")= chr [1:9] "711" "712" "713" "714" ...
#set site order
site_order <- c("S1","S2","S3","S4","S5","S6",
                "S9","S11","S12","S14","S15","S18")  

# Convert Site column to a factor with the custom order
data$Site <- factor(data$Site, levels = site_order)

# Check the structure to confirm the changes
str(data)
## 'data.frame':    1009 obs. of  20 variables:
##  $ Sample.ID            : chr  "3968" "3847" "3841" "3834" ...
##  $ DW.clean.g           : num  30.62 6.59 9.7 5.73 5.73 ...
##  $ DW.wax.g             : num  31.88 6.96 10.22 6.16 6.16 ...
##  $ BW.wax.g             : num  16.16 3.35 5.03 2.89 2.89 ...
##  $ Temperature.C        : num  20.1 20.2 20.2 20.1 20.1 ...
##  $ Notes.x              : chr  "" "" "" "duplicated, Pierrick has page" ...
##  $ Vol.enc              : num  15.72 3.6 5.19 3.27 3.27 ...
##  $ Bulk.Density         : num  1.95 1.83 1.87 1.75 1.75 ...
##  $ Sanger.Plate.Name    : chr  NA "RAA_241009" "RAA_241009" "RAA_241007B" ...
##  $ Sequence.ID          : chr  NA "RA0501" "RA0492" "RA0304" ...
##  $ Plate.Well           : chr  NA "B4" "A3" "C3" ...
##  $ gDNA.Extraction.Plate: chr  NA "Plate015" "Plate015" "Plate013" ...
##  $ Haplotype            : chr  NA "Haplotype 1a" "Haplotype 1a" "Haplotype 1a" ...
##  $ Notes.y              : chr  NA "" "" "" ...
##  $ Species              : Factor w/ 7 levels "P. tuahiniensis",..: NA 3 3 3 3 3 3 2 3 3 ...
##  $ X                    : chr  NA "" "" "" ...
##  $ Stage                : chr  "" "" "" "" ...
##  $ Timing               : chr  "Post_bleach" "Pre_bleach" "Pre_bleach" "Pre_bleach" ...
##  $ Site                 : Factor w/ 12 levels "S1","S2","S3",..: 11 10 10 10 10 10 10 10 10 10 ...
##  $ Transect             : chr  "T2" "T2" "T2" "T2" ...
##  - attr(*, "na.action")= 'omit' Named int [1:9] 711 712 713 714 715 716 717 718 719
##   ..- attr(*, "names")= chr [1:9] "711" "712" "713" "714" ...
bulk.density.plot.raw <- data %>%
  ggplot(aes(x = Site, y = Bulk.Density, color = Species)) +  # Add 'color = Species' to map colors to the Species variable
  geom_point() +  # Plot the data points
  scale_color_manual(values = species_colors) +  # Apply the manual color scale
  facet_wrap(~ Species) +  # Facet by Species
  theme_bw()+
  labs(title="All Raw Data",
    x = "Species",
    y = "Bulk Density g/cm2"
  )  

bulk.density.plot.raw

#remove outliers above 2.93
filtered_data <- data %>%
  filter(Bulk.Density <2.93, na.rm = TRUE) %>%
  filter(Bulk.Density >0, na.rm = TRUE)

bulk.density.plot.outlier.removed <- filtered_data %>%
  ggplot(aes(x = Site, y = Bulk.Density, color = Species)) +  # Add 'color = Species' to map colors to the Species variable
  geom_point() +  # Plot the data points
  scale_color_manual(values = species_colors) +  # Apply the manual color scale
  facet_wrap(~ Species) +  # Facet by Species
  theme_bw()+
  labs(title="Outlier Removed",
    x = "Species",
    y = "Bulk Density g/cm2"
  )  

bulk.density.plot.outlier.removed

bulk.density.plot.outlier.removed.species <- filtered_data %>%
  ggplot(aes(x = Site, y = Bulk.Density, color = Species, group = Species)) +  # Map color to Species
  geom_point(position = position_dodge(width = 0.9), size = 1, alpha = 0.2) +  # Plot the individual data points
  stat_summary(
    fun = "mean",  # Compute the mean
    geom = "point",  # Add the mean as a point
    size = 3,  # Size of the mean points
    shape = 18,  # Diamond shape for mean points
    aes(color = Species),  # Make the mean points the same color as the species
    position = position_dodge(width = 0.9)  # Dodge means to align with the individual points
  ) + 
  geom_line(
    stat = "summary",  # Use the summary statistics (mean) to plot the line
    fun = "mean",  # Line connects the means
    aes(group = Species, color = Species),  # Ensure the line connects the means of each species
    size = 1,  # Thickness of the line
  position = position_dodge(width = 0.9)) +
  scale_color_manual(values = species_colors) +  # Use your manual color palette for the species
  theme_classic() +  # Classic theme for cleaner look
  labs(
    title = "Outlier Removed",
    x = "Site",
    y = "Bulk Density g/cm²"
  )
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Display the plot
#bulk.density.plot.outlier.removed.species

View a single site S1

#remove max outlier
site_filtered_data <- data %>%
  filter(Bulk.Density != max(data$Bulk.Density, na.rm = TRUE)) %>%
  filter(Site == "S1") %>%
  filter(Species != "NA")

bulk.density.plot.S1 <- site_filtered_data %>%
  ggplot(aes(x = Species, y = Bulk.Density, color = Species, group = Species)) +  # Add 'color = Species' to map colors to the Species variable
   geom_point(position = position_dodge(width = 0.9), size = 2, alpha = 0.3) +  # Plot the individual data points
  stat_summary(
    fun = "mean",  # Compute the mean
    geom = "point",  # Add the mean as a point
    size = 5,  # Size of the mean points
    shape = 18,  # Diamond shape for mean points
    aes(color = Species),  # Make the mean points the same color as the species
    position = position_dodge(width = 0.9)  # Dodge means to align with the individual points
  )+
  scale_color_manual(values = species_colors) +  # Apply the manual color scale
  geom_line(
    stat = "summary",  # Use the summary statistics (mean) to plot the line
    fun = "mean",  # Line connects the means
    aes(group = Species, color = Species),  # Ensure the line connects the means of each species
    size = 1,  # Thickness of the line
  position = position_dodge(width = 0.9)) +
  ylim(1.0, 2.4)+
  theme_bw()+
  theme(legend.position ="none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        axis.ticks.x = element_line())+
  labs(title="S1",
    y = "Bulk Density g/cm2"
  )

bulk.density.plot.S1
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

View a single site S2

#remove max outlier
site_filtered_data <- data %>%
  filter(Bulk.Density != max(data$Bulk.Density, na.rm = TRUE)) %>%
  filter(Site == "S2") %>%
  filter(Species != "NA")

bulk.density.plot.S2 <- site_filtered_data %>%
  ggplot(aes(x = Species, y = Bulk.Density, color = Species, group = Species)) +  # Add 'color = Species' to map colors to the Species variable
   geom_point(position = position_dodge(width = 0.9), size = 2, alpha = 0.3) +  # Plot the individual data points
  stat_summary(
    fun = "mean",  # Compute the mean
    geom = "point",  # Add the mean as a point
    size = 5,  # Size of the mean points
    shape = 18,  # Diamond shape for mean points
    aes(color = Species),  # Make the mean points the same color as the species
    position = position_dodge(width = 0.9)  # Dodge means to align with the individual points
  )+
  scale_color_manual(values = species_colors) +  # Apply the manual color scale
  geom_line(
    stat = "summary",  # Use the summary statistics (mean) to plot the line
    fun = "mean",  # Line connects the means
    aes(group = Species, color = Species),  # Ensure the line connects the means of each species
    size = 1,  # Thickness of the line
  position = position_dodge(width = 0.9)) +
  ylim(1.0, 2.4)+
  theme_bw()+
  theme(legend.position ="none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        axis.ticks.x = element_line())+
  labs(title="S2",
    y = "Bulk Density g/cm2"
  )

bulk.density.plot.S2
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

View a single site S3

#remove max outlier
site_filtered_data <- data %>%
  filter(Bulk.Density != max(data$Bulk.Density, na.rm = TRUE)) %>%
  filter(Site == "S3") %>%
  filter(Species != "NA")

bulk.density.plot.S3 <- site_filtered_data %>%
  ggplot(aes(x = Species, y = Bulk.Density, color = Species, group = Species)) +  # Add 'color = Species' to map colors to the Species variable
   geom_point(position = position_dodge(width = 0.9), size = 2, alpha = 0.3) +  # Plot the individual data points
  stat_summary(
    fun = "mean",  # Compute the mean
    geom = "point",  # Add the mean as a point
    size = 5,  # Size of the mean points
    shape = 18,  # Diamond shape for mean points
    aes(color = Species),  # Make the mean points the same color as the species
    position = position_dodge(width = 0.9)  # Dodge means to align with the individual points
  )+
  scale_color_manual(values = species_colors) +  # Apply the manual color scale
  geom_line(
    stat = "summary",  # Use the summary statistics (mean) to plot the line
    fun = "mean",  # Line connects the means
    aes(group = Species, color = Species),  # Ensure the line connects the means of each species
    size = 1,  # Thickness of the line
  position = position_dodge(width = 0.9)) +
  ylim(1.0, 2.4)+
  theme_bw()+
  theme(legend.position ="none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        axis.ticks.x = element_line())+
  labs(title="S3",
    y = "Bulk Density g/cm2"
  )

bulk.density.plot.S3
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

View a single site S4

#remove max outlier
site_filtered_data <- data %>%
  filter(Bulk.Density != max(data$Bulk.Density, na.rm = TRUE)) %>%
  filter(Site == "S4") %>%
  filter(Species != "NA")

bulk.density.plot.S4 <- site_filtered_data %>%
  ggplot(aes(x = Species, y = Bulk.Density, color = Species, group = Species)) +  # Add 'color = Species' to map colors to the Species variable
   geom_point(position = position_dodge(width = 0.9), size = 2, alpha = 0.3) +  # Plot the individual data points
  stat_summary(
    fun = "mean",  # Compute the mean
    geom = "point",  # Add the mean as a point
    size = 5,  # Size of the mean points
    shape = 18,  # Diamond shape for mean points
    aes(color = Species),  # Make the mean points the same color as the species
    position = position_dodge(width = 0.9)  # Dodge means to align with the individual points
  )+
  scale_color_manual(values = species_colors) +  # Apply the manual color scale
  geom_line(
    stat = "summary",  # Use the summary statistics (mean) to plot the line
    fun = "mean",  # Line connects the means
    aes(group = Species, color = Species),  # Ensure the line connects the means of each species
    size = 1,  # Thickness of the line
  position = position_dodge(width = 0.9)) +
  ylim(1.0, 2.4)+
  theme_bw()+
  theme(legend.position ="none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        axis.ticks.x = element_line())+
  labs(title="S4",
    y = "Bulk Density g/cm2"
  )

bulk.density.plot.S4
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

View a single site S5

#remove max outlier
site_filtered_data <- data %>%
  filter(Bulk.Density != max(data$Bulk.Density, na.rm = TRUE)) %>%
  filter(Site == "S5") %>%
  filter(Species != "NA")

bulk.density.plot.S5 <- site_filtered_data %>%
  ggplot(aes(x = Species, y = Bulk.Density, color = Species, group = Species)) +  # Add 'color = Species' to map colors to the Species variable
   geom_point(position = position_dodge(width = 0.9), size = 2, alpha = 0.3) +  # Plot the individual data points
  stat_summary(
    fun = "mean",  # Compute the mean
    geom = "point",  # Add the mean as a point
    size = 5,  # Size of the mean points
    shape = 18,  # Diamond shape for mean points
    aes(color = Species),  # Make the mean points the same color as the species
    position = position_dodge(width = 0.9)  # Dodge means to align with the individual points
  )+
  scale_color_manual(values = species_colors) +  # Apply the manual color scale
  geom_line(
    stat = "summary",  # Use the summary statistics (mean) to plot the line
    fun = "mean",  # Line connects the means
    aes(group = Species, color = Species),  # Ensure the line connects the means of each species
    size = 1,  # Thickness of the line
  position = position_dodge(width = 0.9)) +
  ylim(1.0, 2.4)+
  theme_bw()+
  theme(legend.position ="none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        axis.ticks.x = element_line())+
  labs(title="S5",
    y = "Bulk Density g/cm2"
  )

bulk.density.plot.S5
## Warning: Removed 1 row containing non-finite outside the scale range (`stat_summary()`).
## Removed 1 row containing non-finite outside the scale range (`stat_summary()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

View a single site S6

#remove max outlier
site_filtered_data <- data %>%
  filter(Bulk.Density != max(data$Bulk.Density, na.rm = TRUE)) %>%
  filter(Site == "S6") %>%
  filter(Species != "NA")

bulk.density.plot.S6 <- site_filtered_data %>%
  ggplot(aes(x = Species, y = Bulk.Density, color = Species, group = Species)) +  # Add 'color = Species' to map colors to the Species variable
   geom_point(position = position_dodge(width = 0.9), size = 2, alpha = 0.3) +  # Plot the individual data points
  stat_summary(
    fun = "mean",  # Compute the mean
    geom = "point",  # Add the mean as a point
    size = 5,  # Size of the mean points
    shape = 18,  # Diamond shape for mean points
    aes(color = Species),  # Make the mean points the same color as the species
    position = position_dodge(width = 0.9)  # Dodge means to align with the individual points
  )+
  scale_color_manual(values = species_colors) +  # Apply the manual color scale
  geom_line(
    stat = "summary",  # Use the summary statistics (mean) to plot the line
    fun = "mean",  # Line connects the means
    aes(group = Species, color = Species),  # Ensure the line connects the means of each species
    size = 1,  # Thickness of the line
  position = position_dodge(width = 0.9)) +
  ylim(1.0, 2.4)+
  theme_bw()+
  theme(legend.position ="none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        axis.ticks.x = element_line())+
  labs(title="S6",
    y = "Bulk Density g/cm2"
  )

bulk.density.plot.S6
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

View a single site S9

#remove max outlier
site_filtered_data <- data %>%
  filter(Bulk.Density != max(data$Bulk.Density, na.rm = TRUE)) %>%
  filter(Site == "S9") %>%
  filter(Species != "NA")

bulk.density.plot.S9 <- site_filtered_data %>%
  ggplot(aes(x = Species, y = Bulk.Density, color = Species, group = Species)) +  # Add 'color = Species' to map colors to the Species variable
   geom_point(position = position_dodge(width = 0.9), size = 2, alpha = 0.3) +  # Plot the individual data points
  stat_summary(
    fun = "mean",  # Compute the mean
    geom = "point",  # Add the mean as a point
    size = 5,  # Size of the mean points
    shape = 18,  # Diamond shape for mean points
    aes(color = Species),  # Make the mean points the same color as the species
    position = position_dodge(width = 0.9)  # Dodge means to align with the individual points
  )+
  scale_color_manual(values = species_colors) +  # Apply the manual color scale
  geom_line(
    stat = "summary",  # Use the summary statistics (mean) to plot the line
    fun = "mean",  # Line connects the means
    aes(group = Species, color = Species),  # Ensure the line connects the means of each species
    size = 1,  # Thickness of the line
  position = position_dodge(width = 0.9)) +
  ylim(1.0, 2.4)+
  theme_bw()+
  theme(legend.position ="none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        axis.ticks.x = element_line())+
  labs(title="S9",
    y = "Bulk Density g/cm2"
  )

bulk.density.plot.S9
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

View a single site S11

#remove max outlier
site_filtered_data <- data %>%
  filter(Bulk.Density != max(data$Bulk.Density, na.rm = TRUE)) %>%
  filter(Site == "S11") %>%
  filter(Species != "NA")

bulk.density.plot.S11 <- site_filtered_data %>%
  ggplot(aes(x = Species, y = Bulk.Density, color = Species, group = Species)) +  # Add 'color = Species' to map colors to the Species variable
   geom_point(position = position_dodge(width = 0.9), size = 2, alpha = 0.3) +  # Plot the individual data points
  stat_summary(
    fun = "mean",  # Compute the mean
    geom = "point",  # Add the mean as a point
    size = 5,  # Size of the mean points
    shape = 18,  # Diamond shape for mean points
    aes(color = Species),  # Make the mean points the same color as the species
    position = position_dodge(width = 0.9)  # Dodge means to align with the individual points
  )+
  scale_color_manual(values = species_colors) +  # Apply the manual color scale
  geom_line(
    stat = "summary",  # Use the summary statistics (mean) to plot the line
    fun = "mean",  # Line connects the means
    aes(group = Species, color = Species),  # Ensure the line connects the means of each species
    size = 1,  # Thickness of the line
  position = position_dodge(width = 0.9)) +
  ylim(1.0, 2.4)+
  theme_bw()+
  theme(legend.position ="none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        axis.ticks.x = element_line())+
  labs(title="S11",
    y = "Bulk Density g/cm2"
  )

bulk.density.plot.S11
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

View a single site S12

#remove max outlier
site_filtered_data <- data %>%
  filter(Bulk.Density != max(data$Bulk.Density, na.rm = TRUE)) %>%
  filter(Site == "S12") %>%
  filter(Species != "NA")

bulk.density.plot.S12 <- site_filtered_data %>%
  ggplot(aes(x = Species, y = Bulk.Density, color = Species, group = Species)) +  # Add 'color = Species' to map colors to the Species variable
   geom_point(position = position_dodge(width = 0.9), size = 2, alpha = 0.3) +  # Plot the individual data points
  stat_summary(
    fun = "mean",  # Compute the mean
    geom = "point",  # Add the mean as a point
    size = 5,  # Size of the mean points
    shape = 18,  # Diamond shape for mean points
    aes(color = Species),  # Make the mean points the same color as the species
    position = position_dodge(width = 0.9)  # Dodge means to align with the individual points
  )+
  scale_color_manual(values = species_colors) +  # Apply the manual color scale
  geom_line(
    stat = "summary",  # Use the summary statistics (mean) to plot the line
    fun = "mean",  # Line connects the means
    aes(group = Species, color = Species),  # Ensure the line connects the means of each species
    size = 1,  # Thickness of the line
  position = position_dodge(width = 0.9)) +
  ylim(1.0, 2.4)+
  theme_bw()+
  theme(legend.position ="none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        axis.ticks.x = element_line())+
  labs(title="S12",
    y = "Bulk Density g/cm2"
  )

bulk.density.plot.S12
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

View a single site S14

#remove max outlier
site_filtered_data <- data %>%
  filter(Bulk.Density != max(data$Bulk.Density, na.rm = TRUE)) %>%
  filter(Site == "S14") %>%
  filter(Species != "NA")

bulk.density.plot.S14 <- site_filtered_data %>%
  ggplot(aes(x = Species, y = Bulk.Density, color = Species, group = Species)) +  # Add 'color = Species' to map colors to the Species variable
   geom_point(position = position_dodge(width = 0.9), size = 2, alpha = 0.3) +  # Plot the individual data points
  stat_summary(
    fun = "mean",  # Compute the mean
    geom = "point",  # Add the mean as a point
    size = 5,  # Size of the mean points
    shape = 18,  # Diamond shape for mean points
    aes(color = Species),  # Make the mean points the same color as the species
    position = position_dodge(width = 0.9)  # Dodge means to align with the individual points
  )+
  scale_color_manual(values = species_colors) +  # Apply the manual color scale
  geom_line(
    stat = "summary",  # Use the summary statistics (mean) to plot the line
    fun = "mean",  # Line connects the means
    aes(group = Species, color = Species),  # Ensure the line connects the means of each species
    size = 1,  # Thickness of the line
  position = position_dodge(width = 0.9)) +
  ylim(1.0, 2.4)+
  theme_bw()+
  theme(legend.position ="none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        axis.ticks.x = element_line())+
  labs(title="S14",
    y = "Bulk Density g/cm2"
  )

bulk.density.plot.S14
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

View a single site S15

#remove max outlier
site_filtered_data <- data %>%
  filter(Bulk.Density != max(data$Bulk.Density, na.rm = TRUE)) %>%
  filter(Site == "S15") %>%
  filter(Species != "NA")

bulk.density.plot.S15 <- site_filtered_data %>%
  ggplot(aes(x = Species, y = Bulk.Density, color = Species, group = Species)) +  # Add 'color = Species' to map colors to the Species variable
   geom_point(position = position_dodge(width = 0.9), size = 2, alpha = 0.3) +  # Plot the individual data points
  stat_summary(
    fun = "mean",  # Compute the mean
    geom = "point",  # Add the mean as a point
    size = 5,  # Size of the mean points
    shape = 18,  # Diamond shape for mean points
    aes(color = Species),  # Make the mean points the same color as the species
    position = position_dodge(width = 0.9)  # Dodge means to align with the individual points
  )+
  scale_color_manual(values = species_colors) +  # Apply the manual color scale
  geom_line(
    stat = "summary",  # Use the summary statistics (mean) to plot the line
    fun = "mean",  # Line connects the means
    aes(group = Species, color = Species),  # Ensure the line connects the means of each species
    size = 1,  # Thickness of the line
  position = position_dodge(width = 0.9)) +
  ylim(1.0, 2.4)+
  theme_bw()+
  theme(legend.position ="none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        axis.ticks.x = element_line())+
  labs(title="S15",
    y = "Bulk Density g/cm2"
  )

bulk.density.plot.S15
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

View a single site S18

#remove max outlier
site_filtered_data <- data %>%
  filter(Bulk.Density != max(data$Bulk.Density, na.rm = TRUE)) %>%
  filter(Site == "S18") %>%
  filter(Species != "NA")

bulk.density.plot.S18 <- site_filtered_data %>%
  ggplot(aes(x = Species, y = Bulk.Density, color = Species, group = Species)) +  # Add 'color = Species' to map colors to the Species variable
   geom_point(position = position_dodge(width = 0.9), size = 2, alpha = 0.3) +  # Plot the individual data points
  stat_summary(
    fun = "mean",  # Compute the mean
    geom = "point",  # Add the mean as a point
    size = 5,  # Size of the mean points
    shape = 18,  # Diamond shape for mean points
    aes(color = Species),  # Make the mean points the same color as the species
    position = position_dodge(width = 0.9)  # Dodge means to align with the individual points
  )+
  scale_color_manual(values = species_colors) +  # Apply the manual color scale
  geom_line(
    stat = "summary",  # Use the summary statistics (mean) to plot the line
    fun = "mean",  # Line connects the means
    aes(group = Species, color = Species),  # Ensure the line connects the means of each species
    size = 1,  # Thickness of the line
  position = position_dodge(width = 0.9)) +
  ylim(1.0, 2.4)+
  theme_bw()+
  theme(legend.position ="none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        axis.ticks.x = element_line())+
  labs(title="S18",
    y = "Bulk Density g/cm2"
  )

bulk.density.plot.S18
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

plot to pdf

# Save the plots to a PDF in one column
pdf("Output/BW.plots.pdf", height = 16, width = 10) # Set dimensions for the PDF
grid.arrange(bulk.density.plot.raw,
             bulk.density.plot.outlier.removed,
             bulk.density.plot.outlier.removed.species,
            ncol = 1) # Arrange in 1 column
dev.off() # Close the PDF device
## quartz_off_screen 
##                 2
# Save the plots to a PDF in one column
pdf("Output/BDSite.plots.pdf", height = 18, width = 12) # Set dimensions for the PDF
grid.arrange(bulk.density.plot.S1,bulk.density.plot.S2,
             bulk.density.plot.S3,bulk.density.plot.S4,
             bulk.density.plot.S5,bulk.density.plot.S6,
             bulk.density.plot.S9,bulk.density.plot.S11,
             bulk.density.plot.S12,bulk.density.plot.S14,
             bulk.density.plot.S15, bulk.density.plot.S18,
             ncol = 3) # Arrange in 1 column
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
dev.off() # Close the PDF device
## quartz_off_screen 
##                 2

Make the Moorea Island Map

register_google(key = "AIzaSyAJHWQg-KMSzffFNWaO1zAakoBz-klFhIg") ### use your own API

# location
Moorea<-data.frame(lon = -149.83246425684064, lat = -17.531092816791094)

#Map base
M1<-get_map(Moorea,zoom = 12, maptype = 'satellite')
## ℹ <https://maps.googleapis.com/maps/api/staticmap?center=-17.531093,-149.832464&zoom=12&size=640x640&scale=2&maptype=satellite&language=en-EN&key=xxx-KMSzffFNWaO1zAakoBz-klFhIg>
bbx <- c(left=-149.802,bottom= -17.480,right=-149.805,top=-17.475)
x <- c(bbx["left"], bbx["left"], bbx["right"], bbx["right"])
y <- c(bbx["bottom"], bbx["top"], bbx["top"], bbx["bottom"])
df <- data.frame(x, y)

Mooreamap<-ggmap(M1)+
  scalebar(x.min = -149.90, x.max = -149.05,y.min = -17.4, y.max = -18.0,
           model = 'WGS84', box.fill = c("yellow", "white"), st.color = "white",
           location =  "bottomright", transform = TRUE, dist_unit = "km", dist = 10) +
  #geom_polygon(aes(x=x, y=y), data=df, color="yellow", fill=NA) +
  #geom_text(data = isledata, aes(x=Long, y=Lat, label=island),vjust =0,size=4, color = 'yellow')+
  ggtitle('A')+
  xlab("")+
  ylab("")

Mooreamap
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_text()`).

Make Colony and Transect Map

sitedata <- read.csv("Data/transect.site.coords.csv")
str(sitedata)
## 'data.frame':    12 obs. of  3 variables:
##  $ long     : num  -150 -150 -150 -150 -150 ...
##  $ lat      : num  -17.5 -17.6 -17.6 -17.6 -17.5 ...
##  $ site.name: chr  "S1" "S11" "S12" "S14" ...
labels <- sitedata$site.name

#Map base
#M3<-get_map(NS,zoom = 18, maptype = 'satellite')

MooreaSitemap<-ggmap(M1)+
  scalebar(x.min = -149.75, x.max = -149.79,y.min = -17.60, y.max = -17.58,
           model = 'WGS84', box.fill = c("white", "white"), st.color = "white",st.dist=0.3,
           location =  "bottomright", transform = TRUE, dist_unit = "km", dist = 1) +
  geom_point(data = sitedata, mapping = aes(x=long, y=lat), size=2, color="yellow")+
  geom_text(data = sitedata, aes(x=long, y=lat, label=site.name),vjust = -2,size=3, color="yellow")+
  #geom_line(data = transects, mapping = aes(x=long, y=lat,group=grouping), size=0.5, color="black")+
  #ggtitle('C')+
  xlab("")+
  ylab("")

MooreaSitemap

spatial mapping of plots around the island map

empty_plot <- ggplot() + theme_void()

legend_plot <- ggplot(data = data.frame(x = 1, y = 1, group = names(species_colors)), aes(x = x, y = y, color = group)) +
  geom_point(size = 0) +  # No actual points, just to display the legend
  scale_color_manual(values = species_colors) +  # Use your custom color palette
  theme_void() +  # Remove axes and background
  theme(legend.title = element_blank()) +  # Remove legend title
  guides(color = guide_legend(override.aes = list(size = 5)))
legend_plot

grid_layout <- plot_grid(
  bulk.density.plot.S1,bulk.density.plot.S2,bulk.density.plot.S3,bulk.density.plot.S4,bulk.density.plot.S5,
  bulk.density.plot.S18,empty_plot,empty_plot,empty_plot,bulk.density.plot.S6,
  bulk.density.plot.S15,empty_plot,empty_plot,empty_plot,bulk.density.plot.S9,
  bulk.density.plot.S14, empty_plot,empty_plot,empty_plot,bulk.density.plot.S11,
  empty_plot,empty_plot, bulk.density.plot.S12,empty_plot,legend_plot,
  ncol = 5, nrow = 5
)
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
final_grid_with_map <- grid_layout + 
  draw_plot(MooreaSitemap, x = 0.2, y = 0.2, width = 0.6, height = 0.6)  # Adjust size and position
final_grid_with_map 

ggsave("Output/grid_layout.pdf", final_grid_with_map, width = 10, height = 10)